Trapping of lattice polarons by impurities 
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We consider the effects of single impurities on polarons in three-dimensions (3D) using a continuous 
time quantum Monte-Carlo algorithm. An exact treatment of the phonon degrees of freedom leads 
to a very efficient algorithm and we are able to compute the polaron dynamics on an infinite lattice 
using an auxiliary weighting scheme. The magnitude of the impurity potential, the electron-phonon 
coupling and the phonon frequency are varied. We determine the magnitude of the impurity potential 
required for polaron trapping. For small electron-phonon coupling the number of phonons increases 
dramatically on trapping. The polaron binding diagram is computed, showing that intermediate- 
coupling low-phonon-frequency polarons are localized by exceptionally small impurities. 

PACS numbers: 71.38.-k 



Interest in the role of electron-phonon interactions 
(EPIs) and polaron dynamics has recently gone through 
a vigorous revival. Polarons have been shown to be rele- 
vant in high-temperature superconductors, colossal mag- 
netoresistance oxides, polymer and many inorganic semi- 
conductors, and electron transport through nanowires 
often depends on vibronic displacements of ions. The 
continued interest in polarons extends well beyond the 
physical description of advanced materials. The field has 
been a testing ground for analytical, semi-analytical, and 
different numerical techniques [ij . 

The situation as regards polaron formation and dy- 
namics in real materials is complicated by an intrinsic 
disorder. Moreover, the EPI cannot be considered either 
weak or strong in many of the materials described above, 
so standard approximations based on perturbation the- 
ories break down. In a pioneering paper Economou and 
coauthors ^ studied a one-dimensional (ID) large po- 
laron with diagonal disorder by using methods from the 
theory of non- linear systems. Bronold et al. [3] inves- 
tigated the dynamics of a single electron in a Holstein 
model with a site-diagonal, binary-alloy-type disorder by 
applying a dynamical mean-field theory (DMFT) for a 
Bethe lattice with infinite coordination number. DMFT 
was further improved by Bronold and Fehske 4] , who de- 
scribed Anderson localization and the self-trapping phe- 
nomena within the same model. 

There is a delicate interplay between the self-trapping 
by EPI and trapping by doped-induced disorder in the 
intermediate-coupling regime, so even a single polaron 
has to be studied by exact methods such as continuous- 
time quantum Monte-Carlo (CTQMC) [1,0 and/or dia- 
grammatic Monte-Carlo (DMC) [7!] techniques. In partic- 
ular, the CTQMC algorithm treats the phonon degrees of 
freedom exactly, does not suffer from time discretization 
errors and is not limited to finite lattices, thus providing 
numerically exact solution of the (bi)polaron problem in 
any dimensions and on any lattice. 

In this paper, we use CTQMC to solve the problem 



of the electron interacting with phonons in the presence 
of an impurity on a 3D lattice. The electron hopping 
is assumed to be between nearest neighbors only. The 
phonon subsystem is made up of independent oscillators 
with frequency uj, displacement ^m, momentum Pm = 
—ihd/d£,in and mass M associated with each lattice site. 
The real space Hamiltonian reads. 
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Here (nn') denote pairs of nearest neighbors. The sites 
are indexed by n or m for electrons and ions respec- 
tively. The spin indices and Hubbard U are omitted 
since there is only one electron. The strength of the 
EPI is defined through a dimensionless coupling constant, 
A = /^(0)/2Mw^zt which is the ratio of the polaron 
energy when t = to the kinetic energy of the free elec- 
tron W = zt. In this paper, we discuss Holsteinpolarons, 
which have a force function, /m(n) = K.5n.,m [SJ. A spe- 
cial case of the Hamiltonian given in eq. ([T]) is the case 
where the external potential. An — Af^n.o, representing 
an impurity. We will study this in detail in this paper. 
We note that there is a significant difference between this 
problem, where electrons in both host and impurity are 
coupled to phonons and the traditional electron-phonon 
impurity problem, where only the impurity site is cou- 
pled to a phonon mode [P]. In the problem considered 
here, polarons can exist in both the host and impurity. 

The CTQMC method employed here has been de- 
scribed in detail with regard to the single polaron prob- 
lem in Ref. 5. Here we give a quick overview of the 
extended algorithm for impurities, which is one of the 
main developments here. The effective polaron action 
that results when the phonon degrees of freedom have 
been integrated out analytically is given by. 
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The contribution of this action to the statistical weight 
of a path configuration is e^. $o[r(T), r(r')] = 
Em/m[r(r)]/™[r(T')] (here = uo/t and ^ = t/fcsT) 
is the phonon mediated interaction. 

One of the main compUcations regarding the simula- 
tion of a particle in a single impurity potential is ensuring 
that the whole configuration space is sampled. This can 
be understood in the following way. During a binary kink 
update (discussed below) the ends of the path may either 
stay put or move right or left by one site along one of the 
nearest neighbor bonds, i.e. the path is essentially a ran- 
dom walker. We set up the system with a very small 
attractive A ^ 0~, with the path close to the impurity. 
In ID, the random walker has finite probability to return 
to the start site after a finite number of steps. In 2D the 
walker has vanishing probability to return, so ergodicity 
is not guaranteed. In 3D, the walker will not in general 
return to its start point even after an infinite time, so 
if updates only have the properties of a basic random 
walker, then the Monte-Carlo procedure will fail, since 
the impurity will (on average) never be visited. 

In order to ensure that the path includes sufficient sam- 
pling of the impurity, it would be useful to make the path 
spend more time near the impurity, weighting the esti- 
mators in such a way that the final measurements are 
the same. Sampling in this way may be achieved by 
introducing an auxiliary weighting to the problem. A 
functional of the path w[{C}] in introduced, which rep- 
resents the probability that an unbound path is found 
in a particular configuration, which is tuned so that the 
path samples the impurity regularly. w[{C}] modifies the 
measurements as. 
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and the update probabilities are modified by the ratio of 
weights u;[{D}]/iy[{C}] on changing from configuration 
C to D. Atot is the total action including kinetic energy 
terms (for reasons relating to the continuous time algo- 
rithm, the kinetic energy part of the statistical weight 
is taken into account through the update rule, eq. [6]). 
Thus the measurement of an observable with respect to 
the ensemble defined by the action Atot can be related 



to the measurement of an observable with respect to the 
ensemble defined by Atot and w. Error bar estimation is 
an important aspect of Monte-Carlo simulation. We note 
that it is necessary to take into account the covariance 
between {d/w[{C}]) and {l/w[{C}])j^ „ when 

determining error bars. We use the bootstrap method 
for error analysis. 

The form of w is a matter of choice, but it is useful to 
choose a form that allows control over the confinement 
of the new path. To avoid undesirable long timestep 
correlations, we use an auxiliary weighting of the form 
w oc l/{a + R'^'^^'^^), where R is the distance of the 
configuration from the impurity, d the dimensionality of 
the lattice, 77 is a small value and a stops the weight 
blowing up on the impurity site. In this way, the average 
distance from the impurity is finite, since in the absence 
of interaction, (i?) ~ w[R]R'^~^ dR, shows that "free" 
particles are localized, but not bound too strongly. 

In the path integral QMC for a polaron in an im- 
purity, it is necessary to make update operations in- 
volving two kinks to ensure that the end configura- 
tions of the path remain periodic in imaginary time, 
since the translational symmetry has been broken. A 
binary update satisfying the imaginary-time boundary 
conditions involves adding or removing a kink-antikink 
pair (an antikink to kink 1 is a kink with direction 
—1). Update probability is determined following a sim- 
ilar argument to that in Ref. [iQ] with a small mod- 
ification; consider two path configurations, {C} and 
{D}, where configuration {D} has one more 1 kink at 
time Ti and one more —1 antikink at time T2 than 
{C}. The balance equation is W[{C}]Qa[{C}]P[{C} 
{D}] = W[{D}]Qr[{D}]P[{D} ^ {C}]. With rcla- 
five contribution of configurations {C} and {D} mod- 
ified by the auxifiary weighting, ly [{I?}]/VK[{C}] = 
(tAT)2e^[^^>l-^[^c^>U[{Z?}]/w[{C}]. 

We modify the equal weighting scheme in Ref. [l^ for 
the probability of choosing kinks with a particular kink 
time when adding or removing the second kink, allowing 
weighted kink insertion with the anti-kink time chosen 
with probability p{t ~ r') where p{T)dT = 1. In this 
way, the kinks and antikinks in the pair can be chosen at 
similar r, and the update acceptance is improved. This 
choice leads to update probabilities for insertion of a di- 
rection 1 kink at r and a direction —1 kink at r'. 
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In order to sample the configuration space faster, we 
also examined path updates which shift the whole path 
through a certain distance, but the efficiency of the al- 
gorithm was not improved. The auxiliary weighting ap- 
proach may be used for systems with several particles, 
with the auxiliary weight depending on either the abso- 
lute or relative positions of the paths. 

In the weighted ensemble, the ground state polaron 
energy, number of phonons and average distance from 
the impurity are respectively computed via, 




where N is the total number of kinks in the path. 

Applying the Lang-Firsov canonical transformation 
to the Hamiltonian and assuming large phonon fre- 
quency, an approximate form can be derived, H ~ 
^l^j + J^i^icl^i, with i = texp{-ztX/uj). It 
is possible to compute an analytic value of the impu- 
rity potential, A at which binding occurs. For po- 
larons on a 3D lattice under the influence of an impu- 
rity, Ac = 3.958f = 3.958mo/m*. When there is no 
electron-phonon coupling, the relation Ac = 3.958t is 
exact. While the approximation for t is not expected to 
hold for low phonon frequency, use of the exact value 
of the inverse mass computed from CTQMC in the ab- 
sence of the impurity is expected to lead to a quali- 
tatively correct value for A^. We will return to this 
point later in the article. An exact value for the energy 
of the instantaneous problem is given by the equation, 
1 = I A| / / d'q/{27r)H/[\E\ - 2 ^ti cosfe)]- 

In order to determine the circumstances under which 
the polaron is localized by the impurity, the total energy 
must be determined. Fig. [T] shows the total energy E 
vs A for various A at w = 1. The transition can be seen 
as a sudden change in the gradient. The flat gradient 
corresponds to the energy of the unbound polaron. The 
strong binding asymptotes show how the small polaron 
with large A almost immediately binds strongly as an 
impurity potential in introduced. We show the exact 
energy as a line beneath the A = points, but for the 
other cases, lines are a guide to the eye. 

The measure of the inverse radius gives another crite- 
rion for the critical binding potential. Fig. [2] shows R~^p 
vs A for various A at tj = 1. The inverse radius van- 
ishes continuously on binding. As the electron-phonon 
interaction increases in strength, the polaron binds at 
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FIG. 1: Total energy i5 vs A for various A at a; = 1. In all 
figures, error bars are plotted to 3 standard errors. Where 
error bars are not visible, error bars are smaller than the 
points. The transition can be seen as a sudden change in the 
gradient. The flat gradient corresponds to the energy of the 
unbound polaron. Straight dotted lines indicate the strong 
impurity asymptote —zt\ + A 
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FIG. 2: R^^p vs A for various A at lj = 1. 



monotonically smaller values of the impurity potential, 
until around A = 1, where tiny impurities are capable of 
binding the polaron. This is because the kinetic energy 
of the polaron decreases rapidly on increasing A. Binding 
occurs when the depth of impurity potential is approxi- 
mately the kinetic energy of the free polaron. 

Small (strong-coupling) polarons have much larger 
numbers of phonons associated with the electron than 
large (weak-coupling) polarons. It is therefore of interest 
to see if the reduction in size associated with localization 
is also associated with a similar phenomenon. Fig. [3] 
shows the variation of Nph with A at a number of dif- 
ferent A at oj/i = 1. For small A, there is a sudden 
increase in the number of phonons as the large polaron 
is trapped by the impurity, to small polaron behavior. 
The localization of the polaron wavefunction increases 
the possibility for interaction between the electron and 
the local ion mode through a local coupling. Such an 
effect is expected to be less pronounced for long range 
electron-phonon coupling. 
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FIG. 3: Nph vs A for various \ a,t lj — 1. For small A, there 
is a sudden increase in Nph as the large polaron is trapped by 
the impurity, leading to small polaron behavior. 
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FIG. 4: Binding diagram in A, A space. Also shown are the 
approximate values predicted by the Lang-Firsov transforma- 
tion with zero excited phonons Ac(A) = Ac(0) exp(— zA/o)), 
from the exact mass in the case with translational invariance 
(measured using QMC) Ac (A) — Ac(0)mo/m* and the value 
measured using Monte-Carlo. 

Finally, we show the binding diagram in A, A space 
in Fig. 3) Also shown are the approximate values pre- 
dicted by the Lang-Firsov transformation with zero ex- 
cited phonons Ac(A) = Ac(0) exp(— zA/w), from the ex- 
act mass in the translationally-invariant case (measured 
using QMC) Ac (A) = Ac(0)mo/m* and the value mea- 
sured from the current Monte-Carlo code. For = 1, we 
examined the energy for f3 = 112 close to binding. Below 



f3 = 56, error bars dominated temperature corrections, 
and led to small differences between analytics and nu- 
merics when Lu = 10. At low phonon frequencies, the Hol- 
stein polaron binds almost immediately at A ^ 1. This 
leads to the expectation that polarons with local inter- 
action are almost completely localized in materials with 
strong electron-phonon coupling. Mobile behavior could 
re-emerge if the electron density is similar to the density 
of impurities and a strong Coulomb repulsion {U > 4A) 
is also available so that polarons become repulsive rather 
than attractive impurities when pinned. Longer range 
interactions lead to more mobile polarons, and therefore 
less likely to be pinned (the role of mass is shown by the 
similarity between the binding diagram from QMC, and 
the approximate one determined from the mass of the 
unbound polaron). Such long range electron-phonon in- 
teractions are difficult to justify in 3D materials, and we 
expect that polarons strongly coupled with low-frequency 
phonons are strongly localized in such materials. 

In summary, we have developed an algorithm for study- 
ing the trapping of polarons by impurities. Analysis of 
the total energy showed that for electron-phonon cou- 
pling A > 1, polarons are strongly bound to the impu- 
rity. For small A, there is a more gradual binding, cou- 
pled with a sudden increase in the number of phonons 
present in the polaron. We computed the binding di- 
agram, showing that critical impurity strength changes 
monotonically with A for intermediate-coupling, in con- 
trast to the conclusion for continuum (weak-coupling) po- 
larons, and differs significantly from the strong-coupling 
approximation. 
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